Plotting multispectral data¶
Multispectral data can be plotted as:
- Individual bands
- Spectral indices
- True color 3-band images
- False color 3-band images
Spectral indices and false color images can both be used to enhance images to clearly show things that might be hidden from a true color image, such as vegetation health.
%store -r band_dict ndvi_da den_redlining_gdf
Add missing libraries to the imports
import cartopy.crs as ccrs # CRSs
# Interactive tabular and vector data
import hvplot.xarray # Interactive raster
# Overlay plots
import numpy as np # Adjust images
import xarray as xr # Adjust images
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
import hvplot.pandas
There are many different ways to represent geospatial coordinates, either spherically or on a flat map. These different systems are called Coordinate Reference Systems.
To make interactive geospatial plots, at the moment we need everything to be in the Mercator CRS.
- Reproject your area of interest with
.to_crs(ccrs.Mercator()) - Reproject your NDVI and band raster data using `.rio.reproject(ccrs.Mercator())
Plotting raster and vector data together using pandas
and xarray requires the matplotlib.pyplot
library to access some plot layour tools. Using the code below as a
starting point, you can play around with adding:
- Labels and titles
- Different colors with
cmapandedgecolor - Different line thickness with
line_width
See if you can also figure out what vmin,
robust, and the .set() methods do.
import xarray as xr
#repjecting
redlining_plot_gdf = den_redlining_gdf.to_crs(ccrs.Mercator())
#reproject raster data with raster in and out
ndvi_plot_da = ndvi_da.rio.reproject(ccrs.Mercator())
#reproj band dict
band_plot_dict = {
name: da.rio.reproject(ccrs.Mercator())
for name, da in band_dict.items()
}
ndvi_plot_da.plot(vmin=0, robust=True)
redlining_plot_gdf.plot(ax=plt.gca(), color='none')
plt.gca().set(
xlabel='', ylabel='', xticks=[], yticks=[]
)
plt.title("NDVI showing vegetation health in Denver", fontsize=16, color='green', loc='center')
#no axis values
plt.xlabel("none", fontsize=12, color="blue", labelpad=10) # x-axis label
plt.ylabel("none", fontsize=12, color="green", labelpad=10) # y-axis label
plt.legend()
#desc
# plt.text(2, 20, "NDVI values dark is healthy vegetation and green to yellow is unhealthy to no vegetation.", fontsize=10, color="blue")
plt.show()
/tmp/ipykernel_1626/2217805051.py:12: UserWarning: No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument. plt.legend() /opt/conda/lib/python3.11/site-packages/IPython/core/pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data. fig.canvas.print_figure(bytes_io, **kw)
vegetation along river shows 0.0 - 0.1 healthy vegeation, indication of healthy vegetation along stream. The yellow and green is within city indicating little to no vegetation within the city of denver due to building and other structures
Now, do the same with hvplot. Note that some parameter
names are the same and some are different. Do you notice any physical
lines in the NDVI data that line up with the redlining boundaries?
(
ndvi_plot_da.hvplot(
geo=True,
xaxis=None, yaxis=None
)
* redlining_plot_gdf.hvplot(
geo=True, crs=ccrs.Mercator(),
fill_color=None)
)
/opt/conda/lib/python3.11/site-packages/dask/dataframe/__init__.py:42: FutureWarning: Dask dataframe query planning is disabled because dask-expr is not installed. You can install it with `pip install dask[dataframe]` or `conda install dask`. This will raise in a future version. warnings.warn(msg, FutureWarning)
The following code will make a three panel plot with Red, NIR, and Green bands. Why do you think we aren’t using the green band to look at vegetation?
raster_kwargs = dict(
geo=True, robust=True,
xaxis=None, yaxis=None
)
(
(
band_plot_dict['red'].hvplot(
cmap='Reds', title='Red Reflectance', **raster_kwargs)
+ band_plot_dict['nir'].hvplot(
cmap='Greys', title='NIR Reflectance', **raster_kwargs)
+ band_plot_dict['green'].hvplot(
cmap='Greens', title='Green Reflectance', **raster_kwargs)
)
* redlining_plot_gdf.hvplot(
geo=True, crs=ccrs.Mercator(),
fill_color=None)
)
Red is showing red reflectance which represents unhealthy vegetation since healthy vegetation absorbs redlight on the spectrum, and thus will not reflect red. The red indicates unhealthy to no vegeation. The NIR Reflectance should show healthy vegeation where NIR is high reflectance. Along the river should have high NIR reflectance, however, shows low reflectance and not sure why. This also is the case with the green and red refelectance. Along the river should show low red reflectance but in this image shows high red reflectance. Green indicates healthy vegetation since it is being reflected back in thisimage. You can see the healthy vegetation in along the river.
The following code will plot an RGB image using both matplotlib and hvplot. It also performs an action called “Contrast stretching”, and brightens the image.
- Read through the
stretch_rgbfunction, and fill out the docstring with the rest of the parameters and your own descriptions. You can ask ChatGPT or another LLM to help you read the code if needed! Please use the numpy style of docstrings - Adjust the
low,high, andbrightennumbers until you are satisfied with the image. You can also ask ChatGPT to help you figure out what adjustments to make by describing or uploading an image.
rgb_da = (
xr.concat(
[
band_plot_dict['red'],
band_plot_dict['green'],
band_plot_dict['blue']
],
dim='rgb')
)
def stretch_rgb(rgb_da, low, high, brighten):
"""
Short description:
True color Satelite image of the city of Denver.
Long description...
Parameters
----------
rgb_da: array-like
Where the values are set to show contrast low (1), high(99) and brightest(.2)
param2: ...
...
Returns
-------
rgb_da: array-like
return rgb_da sends the value to produce the output satelite real color image below.
"""
p_low, p_high = np.nanpercentile(rgb_da, (low, high))
rgb_da = (rgb_da - p_low) / (p_high - p_low) + brighten
rgb_da = rgb_da.clip(0, 1)
return rgb_da
rgb_da = stretch_rgb(rgb_da, 1, 99, .2)
rgb_da.plot.imshow(rgb='rgb')
rgb_da.hvplot.rgb(geo=True, x='x', y='y', bands='rgb')
help(stretch_rgb)
Help on function stretch_rgb in module __main__:
stretch_rgb(rgb_da, low, high, brighten)
Short description:
True color Satelite image of the city of Denver.
Long description...
Parameters
----------
rgb_da: array-like
Where the values are set to show contrast low (1), high(99) and brightest(.2)
param2: ...
...
Returns
-------
rgb_da: array-like
return rgb_da sends the value to produce the output satelite real color image below.
Now, plot a false color RGB image. CIR images have the following bands:
- red becomes NIR
- green becomes red
- blue becomes green
You may notice that the NIR band in this image is very bright. Can you adjust it so it is balanced more effectively by the other bands?
rgb_da = (
xr.concat(
[
band_plot_dict['nir'],
band_plot_dict['red'],
band_plot_dict['green']
],
dim='rgb')
)
def stretch_rgb(rgb_da, low, high, brighten):
"""
Short description:
false color Satelite image of the city of Denver.
Long description...
Parameters
----------
rgb_da: array-like
Where the values are set to show contrast low (1), high(99) and brightest(.2)
param2: ...
...
Returns
-------
rgb_da: array-like
return rgb_da sends the value to produce the output satelite false color image below.
"""
p_low, p_high = np.nanpercentile(rgb_da, (low, high))
rgb_da = (rgb_da - p_low) / (p_high - p_low) + brighten
rgb_da = rgb_da.clip(0, 1)
return rgb_da
rgb_da = stretch_rgb(rgb_da, 1, 99, .2)
rgb_da.plot.imshow(rgb='rgb')
rgb_da.hvplot.rgb(geo=True, x='x', y='y', bands='rgb')
help(stretch_rgb)
Help on function stretch_rgb in module __main__:
stretch_rgb(rgb_da, low, high, brighten)
Short description:
false color Satelite image of the city of Denver.
Long description...
Parameters
----------
rgb_da: array-like
Where the values are set to show contrast low (1), high(99) and brightest(.2)
param2: ...
...
Returns
-------
rgb_da: array-like
return rgb_da sends the value to produce the output satelite false color image below.
%store den_redlining_gdf ndvi_da
Stored 'den_redlining_gdf' (GeoDataFrame) Stored 'ndvi_da' (DataArray)
%%capture
%%bash
jupyter nbconvert redlining-42-tree-model.ipynb --to html
--------------------------------------------------------------------------- CalledProcessError Traceback (most recent call last) Cell In[12], line 1 ----> 1 get_ipython().run_cell_magic('bash', '', 'jupyter nbconvert redlining-34-spectral-plot.ipynb --to html\n') File /opt/conda/lib/python3.11/site-packages/IPython/core/interactiveshell.py:2541, in InteractiveShell.run_cell_magic(self, magic_name, line, cell) 2539 with self.builtin_trap: 2540 args = (magic_arg_s, cell) -> 2541 result = fn(*args, **kwargs) 2543 # The code below prevents the output from being displayed 2544 # when using magics with decorator @output_can_be_silenced 2545 # when the last Python token in the expression is a ';'. 2546 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False): File /opt/conda/lib/python3.11/site-packages/IPython/core/magics/script.py:155, in ScriptMagics._make_script_magic.<locals>.named_script_magic(line, cell) 153 else: 154 line = script --> 155 return self.shebang(line, cell) File /opt/conda/lib/python3.11/site-packages/IPython/core/magics/script.py:315, in ScriptMagics.shebang(self, line, cell) 310 if args.raise_error and p.returncode != 0: 311 # If we get here and p.returncode is still None, we must have 312 # killed it but not yet seen its return code. We don't wait for it, 313 # in case it's stuck in uninterruptible sleep. -9 = SIGKILL 314 rc = p.returncode or -9 --> 315 raise CalledProcessError(rc, cell) CalledProcessError: Command 'b'jupyter nbconvert redlining-34-spectral-plot.ipynb --to html\n'' returned non-zero exit status 255.